#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TArc.h"
#include "TEllipse.h"
#include "TClonesArray.h"
#include "../Objects/ObjectClass.hh"

using namespace std;

void gamma(){

  int NO = 3;

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/gamma_tmp.root");
  if(!fread){ cout<<"No Such File"<<endl; return;}

  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  TClonesArray *Particlea=new TClonesArray("Particle");
  TClonesArray *Neutrona=new TClonesArray("Particle");
  TClonesArray *Argona=new TClonesArray("Particle");
  TClonesArray *Gammaa=new TClonesArray("Particle");
  int nPhoton;
  int nElectron;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("Particleb",&Particlea);
  tree->SetBranchAddress("Neutronb",&Neutrona);
  tree->SetBranchAddress("Argonb",&Argona);
  tree->SetBranchAddress("Gammab",&Gammaa);
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nElectron",&nElectron);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  TGraph *NeutronXY[32];
  TGraph *NeutronXZ[32];
  TGraph *ArgonXY[32];
  TGraph *ArgonXZ[32];
  TGraph *GammaXY[32];
  TGraph *GammaXZ[32];
  TGraph *GammaYZ[32];
  for(int i=0;i<32;i++){
    NeutronXY[i]=new TGraph();;
    NeutronXZ[i]=new TGraph();;
    ArgonXY[i]=new TGraph();;
    ArgonXZ[i]=new TGraph();;
    GammaXY[i]=new TGraph();
    GammaXZ[i]=new TGraph();
    GammaYZ[i]=new TGraph();
  }

  TH1F *hNPhoton[8];
  TH1F *hNElectron[8];
  for(int i=0;i<8;i++){
    hNPhoton[i]=new TH1F(Form("hPhoton%d",i),Form("hPhoton%d",i),1000,0,100000);
    hNElectron[i]=new TH1F(Form("hElectron%d",i),Form("hElectron%d",i),500,0,50000);
  }

  int NNEUTRON,NARGON,NGAMMA;
  int Counter=0;
  int cCOPH=0;
  int cCORA=0;
  int cPHRA=0;
  int cCO=0;
  int cPH=0;
  int cRA=0;
  int cTRANS=0;

  int cFid=0;
  int cFidInteraction=0;

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    Particlea->Clear();
    Neutrona->Clear();
    Argona->Clear();
    Gammaa->Clear();

    tree->GetEntry(iev);

    bool IsCompt    = false;
    bool IsPhot  = false;
    bool IsRayl    = false;
    bool IsFid = false;

    //----- Neutron
    int ineutron = 0;
    int trackid = -1;
    int istep = 0;
    for(int i=0;i<Neutrona->GetEntriesFast();i++){
      Particle *neutron = (Particle*)Neutrona->At(i);
      if(trackid!=neutron->TrackID){
	trackid = neutron->TrackID;
	istep = 0;
	ineutron++;
	if( iev==NO ){
	  NeutronXY[ineutron-1]->SetPoint(istep,neutron->xpre,neutron->ypre);
	  NeutronXZ[ineutron-1]->SetPoint(istep,neutron->xpre,neutron->zpre);
	}
      }
      if( iev==NO ){
	//	cout<<istep<<" kine = "<<neutron->ekin<<endl;
	cout<<neutron->TrackID<<neutron->processname<<" "<<neutron->ekin<<endl;
	NeutronXY[ineutron-1]->SetPoint(istep+1,neutron->xpos,neutron->ypos);
	NeutronXZ[ineutron-1]->SetPoint(istep+1,neutron->xpos,neutron->zpos);

      }
      istep++;
    }
    if( iev==NO ){
      NNEUTRON = ineutron;
      cout<<"nNeutron = "<<NNEUTRON<<endl;
    }

    //----- Argon
    int iargon = 0;
    trackid = -1;
    istep = 0;
    for(int i=0;i<Argona->GetEntriesFast();i++){
      Particle *argon = (Particle*)Argona->At(i);
      if(trackid!=argon->TrackID){
	trackid = argon->TrackID;
	istep = 0;
	iargon++;
	if( iev==NO ){
	  ArgonXY[iargon-1]->SetPoint(istep,argon->xpre,argon->ypre);
	  ArgonXZ[iargon-1]->SetPoint(istep,argon->xpre,argon->zpre);
	}
      }
      if( iev==NO ){
	ArgonXY[iargon-1]->SetPoint(istep+1,argon->xpos,argon->ypos);
	ArgonXZ[iargon-1]->SetPoint(istep+1,argon->xpos,argon->zpos);
      }
      istep++;
    }
    if( iev==NO ){
      NARGON = iargon;
      cout<<"nArgon = "<<NARGON<<endl;
    }

    //----- Gamma
    int igamma = 0;
    trackid=-1;
    istep = 0;
    cFidInteraction = 0;
    for(int i=0;i<Gammaa->GetEntriesFast();i++){
      Particle *gamma = (Particle*)Gammaa->At(i);
      if(trackid!=gamma->TrackID){
	trackid = gamma->TrackID;
	istep = 0;
	igamma++;
	if( iev==NO ){
	  GammaXY[igamma-1]->SetPoint(istep,gamma->xpre,gamma->ypre);
	  GammaXZ[igamma-1]->SetPoint(istep,gamma->xpre,gamma->zpre);
	  GammaYZ[igamma-1]->SetPoint(istep,gamma->ypre,gamma->zpre);
	}
      }
      if( iev==NO ){
	//	cout<<istep<<" edep = "<<gamma->edep<<endl;
	GammaXY[igamma-1]->SetPoint(istep+1,gamma->xpos,gamma->ypos);
	GammaXZ[igamma-1]->SetPoint(istep+1,gamma->xpos,gamma->zpos);
	GammaYZ[igamma-1]->SetPoint(istep+1,gamma->ypos,gamma->zpos);
      }

      TString volname = gamma->volumename;

      if( gamma->processname!="Transportation"&&volname=="pFid"){ 
	IsFid=true;
	cFidInteraction++;
      }

      if     ( gamma->processname=="compt"&&volname=="pFid" )    IsCompt=true;
      else if( gamma->processname=="phot"&&volname=="pFid"  )    IsPhot=true;
      else if( gamma->processname=="Rayl"&&volname=="pFid"  )    IsRayl=true;
      
      if( gamma->processname!="Transportation"&&
	  gamma->processname!="compt"&&
	  gamma->processname!="phot"&&
	  gamma->processname!="Rayl") cout<<"***"<<endl<<"ProcName = "<<gamma->processname<<endl<<"***"<<endl;;

      istep++;
    }
    if( iev==NO ){
      NGAMMA = igamma;
      cout<<"nGamma = "<<NGAMMA<<endl;
      cout<<"nPhoton = "<<nPhoton<<endl;
    }


    if(IsFid) cFid++;

    //-----
    hNPhoton[0]->Fill(nPhoton);
    hNElectron[0]->Fill(nElectron);
    Counter++;
    if( IsCompt&&IsPhot ){
      hNPhoton[2]->Fill(nPhoton);
      hNElectron[2]->Fill(nElectron);
      cCOPH++;
    }
    else if( IsCompt&&IsRayl ){
      hNPhoton[6]->Fill(nPhoton);
      hNElectron[6]->Fill(nElectron);
      cCORA++;
    }
    else if( IsPhot&&IsRayl ){
      hNPhoton[7]->Fill(nPhoton);
      hNElectron[7]->Fill(nElectron);
      cPHRA++;
    }
    else if( IsCompt ){
      hNPhoton[1]->Fill(nPhoton);
      hNElectron[1]->Fill(nElectron);
      cCO++;
    }
    else if( IsPhot ){
      hNPhoton[4]->Fill(nPhoton);
      hNElectron[4]->Fill(nElectron);
      cPH++;
    }
    else if( IsRayl ){
      hNPhoton[5]->Fill(nPhoton);
      hNElectron[5]->Fill(nElectron);
      cRA++;
    }
    else{
      hNPhoton[3]->Fill(nPhoton);
      hNElectron[3]->Fill(nElectron);
      cTRANS++;
    }
    
  }

  cout<<"nFid = "<<cFid<<endl;

  TCanvas *c1=new TCanvas("c1","c1",1200,600);
  c1->SetGrid();
  c1->Divide(2,1);
  c1->cd(1);
  TH1 *frame1=gPad->DrawFrame(-180,-180,180,180);
  frame1->SetXTitle("X (mm)");
  frame1->SetYTitle("Y (mm)");
  frame1->Draw();
  c1->cd(2);
  TH1 *frame2=gPad->DrawFrame(-180,-650,180,-300);
  frame2->SetXTitle("Y (mm)");
  frame2->SetYTitle("Z (mm)");
  frame2->Draw();
  TLegend *leg1=new TLegend(0.6,0.2,0.85,0.35);
  leg1->SetFillStyle(0);

  for(int igamma=0;igamma<NGAMMA;igamma++){
    c1->cd(1);
    GammaXY[igamma]->SetLineColor(4);
    GammaXY[igamma]->SetMarkerColor(4);
    GammaXY[igamma]->SetMarkerStyle(21);
    GammaXY[igamma]->SetMarkerSize(1);
    GammaXY[igamma]->Draw("pl");
    c1->cd(2);
    GammaYZ[igamma]->SetLineColor(4);
    GammaYZ[igamma]->SetMarkerColor(4);
    GammaYZ[igamma]->SetMarkerStyle(21);
    GammaYZ[igamma]->SetMarkerSize(1);
    GammaYZ[igamma]->Draw("pl");
  }
  TLegendEntry *le=leg1->AddEntry(GammaXY[0],"Gamma","lp");
  le->SetTextColor(4);
  //
  for(int ineutron=0;ineutron<NNEUTRON;ineutron++){
    c1->cd(1);
    NeutronXY[ineutron]->SetLineColor(2);
    NeutronXY[ineutron]->SetMarkerColor(2);
    NeutronXY[ineutron]->SetMarkerStyle(21);
    NeutronXY[ineutron]->SetMarkerSize(1);
    NeutronXY[ineutron]->Draw("pl");
    c1->cd(2);
    NeutronXZ[ineutron]->SetLineColor(2);
    NeutronXZ[ineutron]->SetMarkerColor(2);
    NeutronXZ[ineutron]->SetMarkerStyle(21);
    NeutronXZ[ineutron]->SetMarkerSize(1);
    NeutronXZ[ineutron]->Draw("pl");
  }
  le=leg1->AddEntry(NeutronXY[0],"Neutron","lp");
  le->SetTextColor(2);
  //
  for(int iargon=0;iargon<NARGON;iargon++){
    c1->cd(1);
    ArgonXY[iargon]->SetLineColor(3);
    ArgonXY[iargon]->SetMarkerColor(3);
    ArgonXY[iargon]->SetMarkerStyle(21);
    ArgonXY[iargon]->SetMarkerSize(1);
    ArgonXY[iargon]->Draw("pl");
    c1->cd(2);
    ArgonXZ[iargon]->SetLineColor(3);
    ArgonXZ[iargon]->SetMarkerColor(3);
    ArgonXZ[iargon]->SetMarkerStyle(21);
    ArgonXZ[iargon]->SetMarkerSize(1);
    ArgonXZ[iargon]->Draw("pl");
  }
  le=leg1->AddEntry(ArgonXY[0],"Argon","lp");
  le->SetTextColor(3);
  //
  TArc *aInner=new TArc(0.,0.,133.,45.,315.);
  TLine *lInner=new TLine(133./sqrt(2),-133./sqrt(2),133./sqrt(2),133./sqrt(2));
  aInner->SetFillStyle(0);
  aInner->SetNoEdges(1);
  aInner->Draw();
  lInner->Draw();
  TLine *lFid1=new TLine(-133,-360,-133,-550);
  TLine *lFid2=new TLine(-133,-550,  94,-550);
  TLine *lFid3=new TLine(  94,-550,  94,-360);
  TLine *lFid4=new TLine(  94,-360,-133,-360);

  c1->cd(1);
  aInner->Draw();
  lInner->Draw();
  leg1->Draw();
  c1->cd(2);
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  leg1->Draw();

  TCanvas *c2=new TCanvas("c2","c2",1200,600);
  c2->Divide(2,1);
  c2->cd(1)->SetGrid();
  c2->cd(1)->SetLogy();

  TH1 *frame2_1=gPad->DrawFrame(0,0.1,50000,1e4);
  frame2_1->SetXTitle("# of Photon");
  frame2_1->SetYTitle("NEV");
  frame2_1->Draw();
  int color[8]={1,2,3,4,6,7,95,51};
  char name[8][32]={"All","Compt","CO + PH","Transport","Phot","Rayl","CO + RA","PH + RA"};

  TLegend *leg2_1=new TLegend(0.6,0.6,0.9,0.9);
  for(int i=0;i<8;i++){
    hNPhoton[i]->SetLineColor(color[i]);
    hNPhoton[i]->Draw("same");
    le=leg2_1->AddEntry(hNPhoton[i],name[i],"l");
    le->SetTextColor(color[i]);
  }
  leg2_1->SetFillStyle(0);
  leg2_1->Draw();

  c2->cd(2)->SetGrid();
  c2->cd(2)->SetLogy();

  TH1 *frame2_2=gPad->DrawFrame(0,0.1,5000,1e4);
  frame2_2->SetXTitle("# of Electron");
  frame2_2->SetYTitle("NEV");
  frame2_2->Draw();
  TLegend *leg2_2=new TLegend(0.6,0.6,0.9,0.9);
  for(int i=0;i<8;i++){
    hNElectron[i]->SetLineColor(color[i]);
    hNElectron[i]->Draw("same");
    le=leg2_2->AddEntry(hNElectron[i],name[i],"l");
    le->SetTextColor(color[i]);
  }
  leg2_2->SetFillStyle(0);
  leg2_2->Draw();


  /*
  c2->cd(2);
  c2_2->SetGrid();
  c2_2->SetLogx();
  c2_2->SetLogy();
  TH1 *frame2_2=gPad->DrawFrame(10000,0.1,100000,1e4);
  frame2_2->Draw();
  int color[5]={1,2,3,4,6};
  for(int i=0;i<5;i++){
    hNPhoton[4-i]->SetLineColor(color[4-i]);
    if(i==1)
    hNPhoton[4-i]->Draw("same");
  }
  */

  cout<<"  All = "<<Counter<<endl;
  cout<<"   CO = "<<cCO<<endl;
  cout<<"   PH = "<<cPH<<endl;
  cout<<"CO+PH = "<<cCOPH<<endl;
  cout<<"CO+RA = "<<cCORA<<endl;
  cout<<"PH+RA = "<<cPHRA<<endl;
  cout<<"   RA = "<<cRA<<endl;
  cout<<"Trans = "<<cTRANS<<endl;



}

void change2(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  const int nfile = 7;
  int No = 6;

  TFile *fread;
  char fname[nfile][16]={"100keV","500keV","800keV","1MeV","2MeV","5MeV","10MeV"};
  float   XX[nfile]    ={     0.1,     0.5,     0.8,   1.0,   2.0,   5.0,   10.0};

  TGraph *gRatio[9];
  TH1F *hDepositE[9];
  for(int i=0;i<9;i++){
    gRatio[i]    = new TGraph();
    hDepositE[i] = new TH1F(Form("hDepositE%d",i),Form("hDepositE%d",i),525,0,2100);
  }

  for(int ifile=0;ifile<nfile;ifile++){
    fread=new TFile( Form("/scratch/shinichi/OutputFiles/neutron_%s.root",fname[ifile]) );
    if(!fread){ cout<<"No Such File: "<<fname[ifile]<<endl; return;}

    cout<<"---------- "<<fname[ifile]<<endl;

    TTree *tree=(TTree*)fread->Get("tree");
    
    TClonesArray *JPMTa=new TClonesArray("PMT");
    TClonesArray *MPMT1a=new TClonesArray("PMT");
    TClonesArray *MPMT2a=new TClonesArray("PMT");
    TClonesArray *MPMT3a=new TClonesArray("PMT");
    TClonesArray *MPMT4a=new TClonesArray("PMT");
    TClonesArray *Particlea=new TClonesArray("Particle");
    TClonesArray *Neutrona=new TClonesArray("Particle");
    TClonesArray *Argona=new TClonesArray("Particle");
    TClonesArray *Gammaa=new TClonesArray("Particle");
    int nPhoton;
    int nElectron;
    
    tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
    tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
    tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
    tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
    tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
    tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
    tree->SetBranchAddress("JPMTb",&JPMTa);
    tree->SetBranchAddress("MPMT1b",&MPMT1a);
    tree->SetBranchAddress("MPMT2b",&MPMT2a);
    tree->SetBranchAddress("MPMT3b",&MPMT3a);
    tree->SetBranchAddress("MPMT4b",&MPMT4a);
    tree->SetBranchAddress("Particleb",&Particlea);
    tree->SetBranchAddress("Neutronb",&Neutrona);
    tree->SetBranchAddress("Argonb",&Argona);
    tree->SetBranchAddress("Gammab",&Gammaa);
    tree->SetBranchAddress("nPhoton",&nPhoton);
    tree->SetBranchAddress("nElectron",&nElectron);
    
    int nev = tree->GetEntries();
    cout<<"NEV = "<<nev<<endl;
    
    int NNEUTRON,NARGON,NGAMMA;
    int Counter=0;
    int CounterFid=0;
    int cESISCP=0;
    int cESIS=0;
    int cESCP=0;
    int cISCP=0;
    int cES=0;
    int cESsingle=0;
    int cIS=0;
    int cCP=0;
    int cTRANS=0;

    vector<float> deES;
    vector<float> deIS;
    vector<float> deCP;
    vector<float> deTR;

    
    for(int iev=0;iev<nev;iev++){
      JPMTa->Clear();
      MPMT1a->Clear();
      MPMT2a->Clear();
      MPMT3a->Clear();
      MPMT4a->Clear();
      Particlea->Clear();
      Neutrona->Clear();
      Argona->Clear();
      Gammaa->Clear();
      
      tree->GetEntry(iev);
      
      bool  IsFiducial = false;
      bool   IsElastic[2];
      bool IsInElastic[2];
      bool   IsCapture[2];
      for(int i=0;i<2;i++){
	IsElastic[i]   = false;
	IsInElastic[i] = false;
	IsCapture[i]   = false;
      }
      int nES = 0;      
      int nIS = 0;      
      int nCP = 0;      
      int nTR = 0;      

      //----- Neutron
      int ineutron = 0;
      int trackid = -1;
      int istep = 0;

      float kinE = 0.;
      deES.clear();
      deIS.clear();
      deCP.clear();
      deTR.clear();

      for(int i=0;i<Neutrona->GetEntriesFast();i++){
	Particle *neutron = (Particle*)Neutrona->At(i);

	TString procname = neutron->processname;
	float kine = neutron->ekin;

	if(trackid!=neutron->TrackID){
	  trackid = neutron->TrackID;
	  istep = 0;
	  kinE = kine;
	  ineutron++;
	}

	float dE = kinE - kine;

	if( neutron->volumename=="pFid" ){
	  IsFiducial = true;

	  if( dE<0 ) cout<<iev<<" "<<procname<<" "<<dE<<endl;
	  if(ifile==No)	  hDepositE[0]->Fill(1e3*dE);

	  if     ( procname=="HadronElastic" ){    IsElastic[1]=true;deES.push_back(1e3*dE);nES++;}
	  else if( procname=="inelastic"     ){  IsInElastic[1]=true;deIS.push_back(1e3*dE);nIS++;}
	  else if( procname=="nCapture"      ){    IsCapture[1]=true;deCP.push_back(1e3*dE);nCP++;}
	  else if( procname=="Transportation"){                      deTR.push_back(1e3*dE);nTR++;}
	  else if( procname!="Transportation") cout<<"ProcName = "<<procname<<endl;
	}
	/*
	if     ( procname=="HadronElastic" )    IsElastic[0]=true;
	else if( procname=="inelastic"     )  IsInElastic[0]=true;
	else if( procname=="nCapture"      )    IsCapture[0]=true;
	else if( procname!="Transportation") cout<<"ProcName = "<<procname<<endl;
	*/
	kinE = kine;
	istep++;
      }
      
      //----- Argon
      int iargon = 0;
      trackid = -1;
      istep = 0;
      for(int i=0;i<Argona->GetEntriesFast();i++){
	Particle *argon = (Particle*)Argona->At(i);
	if(trackid!=argon->TrackID){
	  trackid = argon->TrackID;
	  istep = 0;
	  iargon++;
	}
	istep++;
      }
      
      //-----
      
      Counter++;
      int whichbit = 1;
      //      if( 1 ){
      if( IsFiducial ){
	CounterFid++;
      }

      if( IsElastic[whichbit]&&IsInElastic[whichbit]&&IsCapture[whichbit] ){
	cESISCP++;
      }
      else if( IsElastic[whichbit]&&IsInElastic[whichbit] ){
	if( ifile==No ){
	  for(int iES=0;iES<nES;iES++)  hDepositE[3]->Fill(deES[iES]);
	  for(int iIS=0;iIS<nIS;iIS++)  hDepositE[3]->Fill(deIS[iIS]);
	}
	cESIS++;
      }
      else if( IsElastic[whichbit]&&IsCapture[whichbit] ){
	if( ifile==No ){
	  for(int iES=0;iES<nES;iES++)  hDepositE[7]->Fill(deES[iES]);
	  for(int iCP=0;iCP<nCP;iCP++)  hDepositE[7]->Fill(deCP[iCP]);
	}
	cESCP++;
      }
      else if( IsInElastic[whichbit]&&IsCapture[whichbit] ){
	if( ifile==No ){
	  for(int iIS=0;iIS<nIS;iIS++)  hDepositE[8]->Fill(deIS[iIS]);
	  for(int iCP=0;iCP<nCP;iCP++)  hDepositE[8]->Fill(deCP[iCP]);
	}
	cISCP++;
      }
      else if( IsElastic[whichbit] ){
	if( ifile==No ){
	  for(int iES=0;iES<nES;iES++)  hDepositE[1]->Fill(deES[iES]);
	  if(nES==1){
	    for(int iES=0;iES<nES;iES++)  hDepositE[2]->Fill(deES[iES]);
	  }
	}
	cES++;
	if(nES==1)  cESsingle++;
      }
      else if( IsInElastic[whichbit] ){
	if( ifile==No ){
	  for(int iIS=0;iIS<nIS;iIS++)  hDepositE[5]->Fill(deIS[iIS]);
	}
	cIS++;
      }
      else if( IsCapture[whichbit] ){
	if( ifile==No ){
	  for(int iCP=0;iCP<nCP;iCP++)  hDepositE[6]->Fill(deCP[iCP]);
	}
	cCP++;
      }
      else{
	if( ifile==No ){
	  for(int iTR=0;iTR<nTR;iTR++) hDepositE[4]->Fill(deTR[iTR]); 
	}
	cTRANS++;
      }
      
      
    }

    cout<<"fid/all = "<<CounterFid<<" "<<Counter<<endl;    
    gRatio[0]->SetPoint(ifile,XX[ifile], CounterFid/(float)CounterFid);
    gRatio[1]->SetPoint(ifile,XX[ifile],        cES/(float)CounterFid);
    gRatio[2]->SetPoint(ifile,XX[ifile],  cESsingle/(float)CounterFid);
    gRatio[3]->SetPoint(ifile,XX[ifile],      cESIS/(float)CounterFid);
    gRatio[4]->SetPoint(ifile,XX[ifile],     cTRANS/(float)CounterFid);
    gRatio[5]->SetPoint(ifile,XX[ifile],        cIS/(float)CounterFid);
    gRatio[6]->SetPoint(ifile,XX[ifile],        cCP/(float)CounterFid);
    gRatio[7]->SetPoint(ifile,XX[ifile],      cESCP/(float)CounterFid);
    gRatio[8]->SetPoint(ifile,XX[ifile],      cISCP/(float)CounterFid);


    if( ifile==No ){
      cout<<"ES+IS+CP   = "<<cESISCP<<endl;
      cout<<"ES+IS      = "<<cESIS<<endl;
      cout<<"ES+CP      = "<<cESCP<<endl;
      cout<<"IS+CP      = "<<cISCP<<endl;
      cout<<"ES         = "<<cES<<endl;
      cout<<"ES(single) = "<<cESsingle<<endl;
      cout<<"IS         = "<<cIS<<endl;
      cout<<"CP         = "<<cCP<<endl;
      cout<<"TRANS      = "<<cTRANS<<endl;
    }

    //---
    tree     ->Delete();
    JPMTa    ->Delete();
    MPMT1a   ->Delete();
    MPMT2a   ->Delete();
    MPMT3a   ->Delete();
    MPMT4a   ->Delete();
    Particlea->Delete();
    Neutrona ->Delete();
    Argona   ->Delete();
    Gammaa   ->Delete();
  }


  //------------ Draw Part
  int color[9]={1,2,2,3,4,6,7,95,51};
  char name[9][32]={"All","Elastic","Single Elastic","ES + IS",
		    "Transport","Inelastic","Capture","ES + CP","IS + CP"};

  //-----
  TCanvas *c1=new TCanvas("c1","c1",800,600);
  c1->SetGrid();
  TH1 *frame1=gPad->DrawFrame(0,0,10,0.9);
  frame1->SetXTitle("Neutron Energy (MeV)");
  frame1->SetYTitle("Ratio( nEvent / All )");

  TLegend *leg1=new TLegend(0.6,0.6,0.9,0.9);

  for(int i=0;i<9;i++){
    gRatio[i]->SetLineColor( color[i] );
    gRatio[i]->SetLineWidth( 2 );
    if( i==2 ) gRatio[i]->SetLineStyle(7);
    gRatio[i]->SetMarkerColor( color[i] );
    gRatio[i]->SetMarkerStyle( 21 );
    gRatio[i]->SetMarkerSize( 1 );
    gRatio[i]->Draw("lp");
    TLegendEntry *le=leg1->AddEntry(gRatio[i],name[i],"pl");
    le->SetTextColor(color[i]);
  }
  leg1->SetFillStyle(0);
  leg1->Draw();

  //-----

  TCanvas *c2=new TCanvas("c2","c2",800,600);
  TH1 *frame2=gPad->DrawFrame(0,1e-1,500,1e4);
  c2->SetLogy();
  c2->SetGrid();

  frame2->SetTitle( Form("Neutron Energy = %s",fname[No]) );
  frame2->SetXTitle("Deposit Energy (keV)");
  frame2->SetYTitle("NEV");

  for(int i=0;i<9;i++){
    hDepositE[i]->SetLineColor( color[i] );
    if(i==2)     hDepositE[i]->SetLineStyle( 7 );
    hDepositE[i]->SetLineWidth( 2 );
    hDepositE[i]->Draw("same");    
  }
  leg1->Draw();  
  
}


